{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "a1a5d379",
   "metadata": {},
   "source": [
    "<p align=\"center\">\n",
    "    <img src=\"https://github.com/GeostatsGuy/GeostatsPy/blob/master/TCG_color_logo.png?raw=true\" width=\"220\" height=\"240\" />\n",
    "\n",
    "</p>\n",
    "\n",
    "## Heterogeneity Metrics Demonstration\n",
    "\n",
    "### Calculating Heterogeneity Metrics to Summarize Subsurface Heterogeneity Tutorial\n",
    "\n",
    "* demonstrate various mathematical operations for statistical expectation with distributions\n",
    "\n",
    "* the same workflow is and data is demonstrated in [Excel Heterogeneity Metrics](https://github.com/GeostatsGuy/ExcelNumericalDemos/blob/master/Heterogeneity_Metrics.xlsx).\n",
    "\n",
    "#### Michael Pyrcz, Associate Professor, University of Texas at Austin \n",
    "\n",
    "##### [Twitter](https://twitter.com/geostatsguy) | [GitHub](https://github.com/GeostatsGuy) | [Website](http://michaelpyrcz.com) | [GoogleScholar](https://scholar.google.com/citations?user=QVZ20eQAAAAJ&hl=en&oi=ao) | [Book](https://www.amazon.com/Geostatistical-Reservoir-Modeling-Michael-Pyrcz/dp/0199731446) | [YouTube](https://www.youtube.com/channel/UCLqEr-xV-ceHdXXXrTId5ig)  | [LinkedIn](https://www.linkedin.com/in/michael-pyrcz-61a648a1) | [GeostatsPy](https://github.com/GeostatsGuy/GeostatsPy)\n",
    "\n",
    "#### Heterogeneity Metrics\n",
    "\n",
    "##### Applications with Measures of Heterogeneity\n",
    "\n",
    "* Measures of heterogeneity are often applied proxy, approximate measures to indicate reservoir production / performance.  \n",
    "* These measures may be applied to compare and rank reservoirs or reservoir model realizations for a single reservoir.  \n",
    "\n",
    "##### Best Practice\n",
    "\n",
    "* None of these metrics are perfect.  \n",
    "\n",
    "* The best result possible from rigorous flow forecasting applied to good, full 3D reservoir models, use the physics when possible!\n",
    "\n",
    "* Integrate all relevant information, at sufficient scale to resolve important features\n",
    "\n",
    "###### Caution\n",
    "\n",
    "* Use of simple heterogeneity measures for ranking reservoir and reservoir models can be dangerous.  \n",
    "\n",
    "* Inaccuracy can result in incorrect rank estimates; therefore, incorrect business decisions.\n",
    "\n",
    "###### Other Measures\n",
    "\n",
    "* We just consider simple, static measures here\n",
    "\n",
    "* I also have a Python demonstration for [Lorenz coefficient](https://github.com/GeostatsGuy/PythonNumericalDemos/blob/master/Lorenz_coefficient_demo.ipynb)\n",
    "\n",
    "* Got your own measure?  You may develop a new metric. Novel methods for quantifying heterogeneity within reservoirs is a currently active area of research.\n",
    "\n",
    "#### Objective \n",
    "\n",
    "Let's demonstate these heterogeneity metrics on a subsurface dataset.\n",
    "\n",
    "#### Getting Started\n",
    "\n",
    "Here's the steps to get setup in Python:\n",
    "\n",
    "* Install Anaconda 3 on your machine (https://www.anaconda.com/download/).  \n",
    "\n",
    "#### Load the Required Libraries\n",
    "\n",
    "The following code loads the required libraries."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "64adca6c",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np                                                # arrays\n",
    "import pandas as pd                                               # dataframes\n",
    "import scipy.stats as stats                                       # statistical functions\n",
    "import matplotlib.pyplot as plt                                   # plotting\n",
    "import math                                                       # square root to calculate standard deviation from variance"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0300e7e9",
   "metadata": {},
   "source": [
    "#### Set the working directory\n",
    "\n",
    "I always like to do this so I don't lose files and to simplify subsequent read and writes (avoid including the full address each time). "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "314135ac",
   "metadata": {},
   "outputs": [],
   "source": [
    "#os.chdir(\"c:/PGE383/Examples\")             # set the working directory"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f46a8eef",
   "metadata": {},
   "source": [
    "#### Loading Tabular Data\n",
    "\n",
    "Here's the command to load our comma delimited data file in to a Pandas' DataFrame object."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "fb7bd372",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>Depth</th>\n",
       "      <th>Porosity</th>\n",
       "      <th>Perm</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>0.25</td>\n",
       "      <td>12.993634</td>\n",
       "      <td>265.528738</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>0.50</td>\n",
       "      <td>13.588011</td>\n",
       "      <td>116.891220</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>0.75</td>\n",
       "      <td>8.962625</td>\n",
       "      <td>136.920016</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>1.00</td>\n",
       "      <td>17.634879</td>\n",
       "      <td>216.668629</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>1.25</td>\n",
       "      <td>9.424404</td>\n",
       "      <td>131.594114</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "   Depth   Porosity        Perm\n",
       "0   0.25  12.993634  265.528738\n",
       "1   0.50  13.588011  116.891220\n",
       "2   0.75   8.962625  136.920016\n",
       "3   1.00  17.634879  216.668629\n",
       "4   1.25   9.424404  131.594114"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "df = pd.read_csv('https://raw.githubusercontent.com/GeostatsGuy/GeoDataSets/master/PorPermSample1.csv')\n",
    "df.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6ad4e47d",
   "metadata": {},
   "source": [
    "#### Feature Engineering\n",
    "\n",
    "We will need to make a new feature, rock quality index, the ratio of permeability divided by porosity."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "23e763c4",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>Depth</th>\n",
       "      <th>Porosity</th>\n",
       "      <th>Perm</th>\n",
       "      <th>PermPor</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>0.25</td>\n",
       "      <td>12.993634</td>\n",
       "      <td>265.528738</td>\n",
       "      <td>20.435295</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>0.50</td>\n",
       "      <td>13.588011</td>\n",
       "      <td>116.891220</td>\n",
       "      <td>8.602526</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>0.75</td>\n",
       "      <td>8.962625</td>\n",
       "      <td>136.920016</td>\n",
       "      <td>15.276776</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>1.00</td>\n",
       "      <td>17.634879</td>\n",
       "      <td>216.668629</td>\n",
       "      <td>12.286369</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>1.25</td>\n",
       "      <td>9.424404</td>\n",
       "      <td>131.594114</td>\n",
       "      <td>13.963124</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "   Depth   Porosity        Perm    PermPor\n",
       "0   0.25  12.993634  265.528738  20.435295\n",
       "1   0.50  13.588011  116.891220   8.602526\n",
       "2   0.75   8.962625  136.920016  15.276776\n",
       "3   1.00  17.634879  216.668629  12.286369\n",
       "4   1.25   9.424404  131.594114  13.963124"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "df['PermPor'] = df['Perm'].values/df['Porosity'].values\n",
    "df.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "52a22fd8",
   "metadata": {},
   "source": [
    "#### Variance of Permeability\n",
    "\n",
    "It is common to use the sample or population variance. I demonstrate the sample mean. \n",
    "\n",
    "\\begin{equation}\n",
    "\\sigma_{X_k}^2 = \\frac{1}{n} \\sum_{i=1}^{n-1} (x_{k,i} - \\overline{x_k})^2\n",
    "\\end{equation}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "id": "ad95297b",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Sample Variance of Permeability: 6544.83 mD^2\n"
     ]
    }
   ],
   "source": [
    "var_perm = np.var(df['Perm'].values,ddof=1)\n",
    "print('Sample Variance of Permeability: ' + str(np.round(var_perm,2)) + ' mD^2')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5372e732",
   "metadata": {},
   "source": [
    "#### Coefficient of Variation of Permeability\n",
    "\n",
    "Another common metric is the coefficient of variation, the standard deviation standardized by the mean.\n",
    "\n",
    "\\begin{equation}\n",
    "C_{v_k} = \\frac{\\sigma_k}{\\overline{k} }\n",
    "\\end{equation}\n",
    "\n",
    "Note, by specifying the ddof arguement of 1, we are using the sample standard deviation in the calculation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "id": "fa8ac2d4",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Coefficeint of Variation of Permeability: 0.502 unitless\n"
     ]
    }
   ],
   "source": [
    "coefvar_perm = stats.variation(df['Perm'].values,ddof=1)\n",
    "print('Coefficeint of Variation of Permeability: ' + str(np.round(coefvar_perm,3)) + ' unitless')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ab82a81d",
   "metadata": {},
   "source": [
    "#### Coefficient of Variation of Permeability / Porosity\n",
    "\n",
    "Also it is common to calculate the coefficient of variation of the rock quality index feature, permeability divided by porosity.\n",
    "\n",
    "\\begin{equation}\n",
    "C_{v_\\frac{k}{\\phi}} = \\frac{\\sigma_{\\frac{k}{\\phi}}}{\\overline{\\frac{k}{\\phi}} }\n",
    "\\end{equation}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "id": "675edec9",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Coefficient of Variation of Permeability / Porosity: 0.369 unitless\n"
     ]
    }
   ],
   "source": [
    "coefvar_permpor = stats.variation(df['PermPor'].values,ddof=1)\n",
    "print('Coefficient of Variation of Permeability / Porosity: ' + str(np.round(coefvar_permpor,3)) + ' unitless')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "798b1fef",
   "metadata": {},
   "source": [
    "#### Dykstra Parsons\n",
    "\n",
    "Now let's calculate the Dykstra-Parsons coefficient. \n",
    "\n",
    "\\begin{equation}\n",
    "DP = \\frac{P50_k - P16_k}{P50_k}\n",
    "\\end{equation}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "id": "d19fbcdf",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Permeability P16: 84.797\n",
      "Permeability P50: 144.33\n",
      "\n",
      "Dykstra-Parsons Coefficient: 0.412 unitless\n"
     ]
    }
   ],
   "source": [
    "P16_perm = np.percentile(df['Perm'].values,16)\n",
    "P50_perm = np.percentile(df['Perm'].values,50)\n",
    "\n",
    "print('Permeability P16: ' + str(np.round(P16_perm,3)))\n",
    "print('Permeability P50: ' + str(np.round(P50_perm,3)))\n",
    "\n",
    "dp = (P50_perm - P16_perm)/P50_perm\n",
    "print('\\nDykstra-Parsons Coefficient: ' + str(np.round(dp,3)) + ' unitless')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f3d6e017",
   "metadata": {},
   "source": [
    "#### Dykstra-Parsons Improved by Fitting a Lognormal Distribution to Permeability\n",
    "\n",
    "We may be able to improve our Dykstra-Parsons coefficient calculation by fitting a lognormal distibution to permeability and then using the parametric lognormal distribution to get a more accurate estimate of the P50 and P16 or permeability.\n",
    "\n",
    "1. calculate the lognormal parameters, \\mu (mu) and \\sigma (sigma)\n",
    "2. use the parametric lognormal CDF inverse to get the P50 and P16"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "id": "455dcefe",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Lognormal distribution parameters are, mu: 4.97, and sigma: 0.47\n"
     ]
    }
   ],
   "source": [
    "mean_perm = np.average(df['Perm'].values)\n",
    "\n",
    "mu = np.log((mean_perm**2)/math.sqrt(var_perm + mean_perm**2))\n",
    "sigma = math.sqrt(np.log(var_perm/(mean_perm**2)+1))\n",
    "\n",
    "print('Lognormal distribution parameters are, mu: ' + str(np.round(mu,2)) + ', and sigma: ' + str(np.round(sigma,2)))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7b36506b",
   "metadata": {},
   "source": [
    "Now we can calculate the P50 and P16 of permeability from the fit lognormal distribution.\n",
    "\n",
    "* Note, we could have been more rigorous with our distribution fit with a ordinary least squares or maximum likelihood approach."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "16789c60",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Permeability P16 from Lognormal Parametric Distribution: 89.754\n",
      "Permeability P50 from Lognormal Parametric Distribution: 143.869\n"
     ]
    }
   ],
   "source": [
    "P16_perm_lognorm = stats.lognorm.ppf(0.16,s = sigma, scale = math.exp(mu))\n",
    "P50_perm_lognorm = stats.lognorm.ppf(0.50,s = sigma, scale = math.exp(mu))\n",
    "\n",
    "print('Permeability P16 from Lognormal Parametric Distribution: ' + str(np.round(P16_perm_lognorm,3)))\n",
    "print('Permeability P50 from Lognormal Parametric Distribution: ' + str(np.round(P50_perm_lognorm,3)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "762a1c08",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "Dykstra-Parsons Coefficient: 0.376 unitless\n"
     ]
    }
   ],
   "source": [
    "dp_lognorm = (P50_perm_lognorm - P16_perm_lognorm)/P50_perm_lognorm\n",
    "print('\\nDykstra-Parsons Coefficient: ' + str(np.round(dp_lognorm,3)) + ' unitless')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "af1c596d",
   "metadata": {},
   "source": [
    "#### Visualize the Lognormal Permeability Parametric Distribution\n",
    "\n",
    "Let's check our lognormal parametric distribution fit.\n",
    "\n",
    "* we visualize the lognormal parametric CDF and the permeability data nonparametric distribution together to check the fit"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "id": "8c05f34f",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAewAAAGWCAYAAAC6k7DHAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAABYQElEQVR4nO3deXxU1f3/8deHsARkExABkV0sCAouWJVWUGvVKoorS1UUpPxUbGtLxcbWWI1atailKqVGQU2CqEDxq9Yqkla01qq4AHEJAhGRRVAWWUPO749zk0ySyWSATGYm834+Hnlk5t47dz5z5s585p5z7jnmnENEREQSW4N4ByAiIiI1U8IWERFJAkrYIiIiSUAJW0REJAkoYYuIiCQBJWwREZEkoIRdi8xsjJktirD+JTO7Mty2ZrbNzHrURZyS+EKPlTjGMMPM7ohnDInCzIaY2ep4xxEPZvZbM3s0GfYd+j1a28evmU0zs9/V1v72R8ImbDNbaWY7gjdgnZk9bmbN4x3XgXDOne2cm1nNuubOuc/hwA+0kLLbambfmtmbZjbBzKJ6v82sm5k5M2t4ADE4M/so9DnN7A4zm7G/+6wPzCzTzJ6qabtIx0oiqOnHaaoJjvfvgu+rL81sipmlxTuuSKL9EeKcu9M5N24/9p9vZjuD76EtZvaumU02syb7uu9gXzVuF/o9eiDCHd/OuQnOudsPdN8HImETduA851xz4FjgBOCWfXmweYn+GmPlPOdcC6ArcDdwE5BdxzF0AkbU8XPukwP5URILKX7M1po4va/HBN9XpwOjgGv25cGJdixCrcR0ffA91BH4Ff774EUzswMOLkQill0sJMUXg3PuS+AloB+AmX0/OGv81sw+MLMhpdsGv8SyzOwNYDvQI/j1e62ZfRb82rvdzHqa2X+CX36zzaxxyD7ONbP3Q85Ojw5ZN9nMlgf7WWZmwyuFa2Y21cw2m9nHZnZ6pdjC/koMYuxlZuOB0cBvgl/rz5vZJDN7rtL2U83sgSjKbrNzbj5wGXClmZWW4U/MbHHw+r8ws8yQh/07+P9tEMNJQXm9ZmYbzexrM8sxs9Y1PP09wG3VfZjMbJiZLQ3KOd/M+oSsW2lmvzazD4OyfNrM0oN1Q8xstfnqtK+DbUeHPLba12bltQdjzawIeC1Y/oyZrQ2e699mdlTIY2aY2cPmq6m3mdkbZtbBzB4ws2+C93lgyPadzOw5M9tgZivM7IZg+VnAb4HLgv18ECwPd8xWOFbM7BozKwg57o6tpkwfDF5z6RnND0LWZQbH+hPBfpaa2fEh6wea2XvBuqeB9Mhvb3hmdrKZ/S8oy/+Z2ckh67oH5bvVzF41s4csqHEIeW+uNLOi4L3NCHlsk6DM1wR/D1hwthZyTNxkZmuBx4PX+4yZPRU830dm1tvMbjaz9UE5nRmy/6tCyvhzM/vZ/rx+59zHwOtAP6vhcxMcuzeZ2YfAd2bW0CJ8x5g/83vDzO43/7n5PCjvMcHrWW8hTSlBmd0XlOc689W6Tc3sIPx3aqfgWNwWHLeZZvZsUGZbgDFWqVbIzAZb+ffvF2Y2Jooy+c45lw8MA04CfhLsq2zfZpYePO/GYN//M7NDzSwL+AHwlyDOvwTbOzO7zsw+Az4LWdYr5KnbmdkrQVn+y8y6BttVqUW04DNn/ntoGnBS8HzfBusr1Hya/0wWmtkmM5tvZp1C1jnztZqfmf+OeMjM/0gx/z3/L/Ofj6/Nf9ai45xLyD9gJXBGcPtwYClwO3AYsBE4B/+D40fB/UOCbfOBIuAooCHQCHDAfKBlsHwXsADoAbQClgFXBo8/FlgPnAikAVcGsTQJ1l+CP3NsgE+C3wEdg3VjgGLgl8HzXgZsBtqExDYuZNtFIa/XAb2C2zOAO0LWdQyep3Vwv2EQ43E1lV2l5UXA/wtuDwH6B6/jaGAdcEGwrlsQT8OQx/YKyroJcAg+qT8Q4f1zwBHAuyGv+Q5gRnC7d/CafhSU1W+AQqBxyGt4OyjrNkABMCEk9mJgShDPqcG+jtyH1/YEcBDQNFh+NdAi2N8DwPshr2UG8DVwHD6JvQasAK7AHyN3AAuDbRsEr/n3QGP8MfY58ONgfSbwVKWyyqfqMZsfUm6XAF/ia5kseC+6VlPuPwXaBvv5FbAWSA957p34z04acBfwVrCuMbCK8mP3YmAPIcdhpecZQ8jxG7K8DfANcHkQw8jgfttg/X+A+4LnGwxsKS2PkPfmb0BT4Bj8Z7VPsP4PwFtAe/wx+CZwe6Vj4o/Be9g05PX+OIjlieB9ywhe4zXAipDYfwL0DMr4VPyPp2ND9r+6huO99PPbNyj3sdTwucEf5+/jv+NKj8VovmOuovzYKwIeCp7jTGAr0DzY/gH8d18b/PH9PHBXda8pKLM9wAXB85eWY+l71CXY/8igDNsCA6opk3yCY7jS8n8Df6z8eQB+FsTXLHhtxwEtq9tXUOavBK+taZj3YUYQ6w+DsnmQ4Jgl/Hdc2XMQ5vgm5HsZOA3/nXBssO+pwL8rxfZ/QOugzDYAZwXr8vDHYAP898ngaHKicy7hE/Y24Fv8F8nDwcFzE/BkpW1fpjzh5gN/CPPGnhJy/13gppD7fyL4EAGPEHwJhKz/BDi1mjjfB84PeZPXABay/m3g8poOiDAH2h2Vnucl4Jrg9rnAshrKLlzCfgvIqOYxDwD3V3cwh9n+AmBxhPUO/2V1Dv4LpQkVE/bvgNkh2zfAJ6UhIa/hpyHr7wGmBbeH4L+0DgpZPxv43T68th4RYm8dbNMq5P34W8j6iUBByP3+wLfB7ROBokr7uxl4PLidSfiEXfmYDT1WXgZ+vp+fo2/wVbWlz/1qyLq+wI7g9g+peuy+Wfk4DFk3hvAJ+3Lg7UrL/hNs3yV435qFrHuKqgm7c6XPz4jg9nLgnJB1PwZWhhwTuwl+nIS83ldC7p+H/05JC+63CJ6vdTWvcV5puRNdwt4SlPdy/LHeoKbPDf44v7qG9/B9Kn7HfFbp2HPAoSHLNgID8D88vgN6hqw7ieBHSrjXFJTZv8MsK32PbgbmRnns5RM+Yc8i+DxV2vfVwTF3dDT7Cl73aWGWhX6PzgpZ1xzYi/9xVHqs7W/CzgbuqbTvPUC3kDgGh6yfDUwObj8BTCfkOI/2L9GrxC9wzrV2znV1zl3rnNuBb5O9JKgy+TaorhiMPwst9UWYfa0Lub0jzP3SDm1dgV9V2v/h+F+8mNkVVl5d/i2+mr5dyL6+dMG7ElhV+tgDNBN/9kTw/8n92MdhwCYAMzvRzBaar7bdDEyg4uuowMzam9ks8x1qtuC/aKvdvpRz7kV8wh5faVUnfNmUbleCf98OC9lmbcjt7ZS/RwDfOOe+C7lfVs5RvrayY8TM0szs7qAacgv+S5RKj9mX46dTpePnt8ChRBbumC11OD4J1MjMfmW+Wndz8NytqPg6KpdpelAt2Inwx+6+qvC+huznsGDdJufc9pB14V53de975X1X/mxtcM7trLSvyu/T1865vSH3Kd2/mZ1tZm8FVZzf4n9s1niMhzjWOXewc66nc+4W51xJlJ+bCmUQxXdM5deEcy7c8XgI/mz13ZB9/SNYHkmtHIsRlH0PVfIk/sfpLPNNHveYWaMa9hUp1grrnXPbguetje/jyt9f2/A/lKL5/voN/sfU2+abpa6O9kkTPWGH8wX+DLt1yN9Bzrm7Q7Zx1T04yv1nVdp/M+dcXtD+8TfgenwVX2tgCb7wSx1W2lYR6II/c9kX4eKfBxxtvg36XCBnX3ZoZifgD6bSno+5+Kqyw51zrfBtNqVxh3v+u4LlRzvnWuJ/NETbceQWfBVQs5Bla/DJrTQ+w38ZfBnlPg8O2uFKhZZzpNdWKvQ1jgLOB87AJ7hupWFFGUuoL/BnMKHHTwvn3Dlhnre6eMLts2dNT2y+vfom4FLg4OD43Ex0r+Mrwh+7+6rC+xqyny+D52hjZqHHweEHsO/Kn639/tybbwt/Dl9df2hQdi+yf8dAqGg+N2VxR/kdE62v8cn7qJBjsZXzHeMqPG918YQR1bFYHTM7HF/V/XqVJ3Vuj3PuNudcX+Bk/PfcFQcQK4QcX+avMmqDP2ZKf+yHHosd9mG/lb+/DsI3D9T4/eWcW+ucu8Y51wnfDPBwpXb3aiVjwn4KOM/MfhycGaWb73DSuZb2/zdgQnCWZmZ2kPlOTC3wbZ4O3x6BmV1F0BEuRHvgBjNrZGaXAH3wH/x9sQ7f9lkmOHN4Fp+M3nbOFUWzIzNraWbn4quhnnLOfRSsaoE/29lpZoPwSavUBqCkUgwtCJoozOwwYFK0L8b5ziYf4fsDlJoN/MTMTg9+Rf8K3175ZrT7xXdoaxwkqnOBZ6J4beG0CJ57I/4DfOc+xFDZ28AW8x2JmgbHaL/gBxP497ab7VtP8EeBX5vZccEx2Sv4Yg/3Oorx719DM/s9vt9GNP4TPPYG8x2fLgQG1fAYCz5/ZX/4Y723mY0K9nMZvur9/5xzq4B3gMzgfTsJX00drTzgFjM7xMza4fsJ1HiJXJQa45ttNgDFZnY2vj34QO3r5yaa75ioBLVWfwPuN7P2wf4OM7MfB5usA9qaWat92G0OcIaZXRq8v23NbEBNDzKzZmZ2KvB3/GekyneimQ01s/7mL4fbgq9iLq0NqfKdGKVzzHeSa4zvA/Vf59wXzrkN+OT60+AzejUVf4isAzpbSGfkSnKBq8xsQPBj785g3ytrCsjMLgnJV9/g3++9ER5SJukStnPuC/zZ0G/xB/UX+A9BrbwW59w7+M4of8EXZiG+PQPn3DJ8e/d/8G9of+CNSrv4L76z1ddAFnCxc27jPoaRDfQNqrHmhSyfGTxnNNXhz5vZVnz5ZOA7aF0Vsv5a4A/BNr/HJ1AAgirLLOCNIIbvA7fhO1hsBl4A5uzja7oF/+u29Dk+wZ9tTMWX1Xn4S9F2R7m/tfj3Zw3+S2SC871zI762ajyBr976Et8B8a0oY6giqG49D9+GuAL/2h7Fn7lD+Y+KjWb2XpT7fAb/fuTiO9HMI6QsQ7yM7+vwKf717KTmKsPS59gNXIg/1r/Bd3aq6T0+GX8GF/q3Gf/j6Vf4H0C/Ac51zn0dPGY0vh11I76d92n8j6Vo3IFP+B/ifwC+Fyw7YM65rcAN+GPlG/yPvPm1sOt9+txE+R2zL27Cf4e9FVTJvwocGTzXx/gfQZ8Hn/Maq4qDE4Vz8O/vJnz7+jERHvKX4HO4Dt+X5Dl856uSMNt2wJ+UbMF3Mv0X5T/IHgQuNt/j+s81xRkiF7g1iPU4/PFX6hp87tiI7/AZerLwGr6j81oz+5pKnHML8P1wnsPXHPUk+ktYTwD+a2bb8MfYz51zK6J5oFVsspJEZmZdgI+BDs65LfGOJ17MX8b3lHOutmpVJE7MX9LysXPu1njHIpLoku4MO1UFVag34ns9pmyyluRmZieYvza5gfnr0s/H1xiISA1ilrDN7DHzF/EvqWa9mdmfzV94/qFVMxCElHVo2IK/nlNnIpLMOuAvn9kG/Bk/LsDiuEYkkiRiViVuZj/EfyifcM5V6TRhZufgr2c9B3/t6oPOuRNjEoyIiEiSi9kZtnPu34S/1q7U+fhk7pxzbwGtzaxjhO1FRERSVjwHTD+Mij1YVwfLvqq8ofnxtccDpKenH9ely/5cHpoaSkpKaNBAXROqo/KJTOUTmconsqQqn8q1y5Xub92yhZ3fbeLb74ppfVBD0g86mBbNW0T3+GrWfVpU9LVzrqaBa6oVz4QdbiCAsPXzzrnp+KHcOPLII90nn3wSy7iSWn5+PkOGDIl3GAlL5ROZyicylU81nIM9e8hftIghxxwDu3f7v127ym+Hux/NNvt6v7pt9uyJ+uXkARnNIedyGNwdFq0oZuxTG8hcv4GRoRuaQZMm0Lix/yu9Xfl/8GdFRfszemCZeCbs1VQc5agz+z4imIhI/ebcvie5ukyEu3eXJ8P77oPTT4/8eqJVmgwrJb1wiZDWrSOvj2YfIfezrhpF9jnrGBqMPza0F2T/FCa+0pmR73xU/piG+5hCD3BW0Xgm7PnA9WY2C9/pbLNzrkp1uIhIUiopgQ0b4Kuvwv99/XV5wouULPfhzDBqDRrUnNRK7zdrFn3S69wZ/vzn2kmcaWm1/7qjVLB6PYO7V1w2uDsUrPjS/ziIk5glbDPLw88G087MVuMvR2oE4Jybhh+a7hz8KDzbqTgKl4hIrcvLzSErM4OC5UX06dmFjMwsRo4aXfMDQ+3ZA+vWVUy+a9ZUTcjr1sHeMCNOHnwwdOwIhxziv/yjTZz7sk1N92OVDPPz4bLLYrPvOtSnZxcWrVhVdoYNsGiFXx5PMUvYzrmRNax3wHW18Vx79uxh9erV7NxZeaKe1NOqVSsKCgriHUbC2tfySU9Pp3PnzjRqVNOkQZLo8nJzyLhxPNnDtzN4HCxasYqxN44H8El7586ak/DVV8Npp1XtVGTmE3DHjv7v6KP9/06dypd17AgdOkB6ehxeveyLjMwsxpYeK919sh47txlZU7LiGlc8q8RrzerVq2nRogXdunXDDrCNINlt3bqVFi1a1LxhitqX8nHOsXHjRlavXk337t1rfoAkrq1byfrtr8kevr1iu+Tw7UwcP4aR110P335b9XFpaT7JduwIXbv6s+Nbb62YhDt2hPbtQT/q6o3SWpeJIbUxWVP2ozamltWLhL1z504la6l1Zkbbtm3ZsGFDvENJaLVSzVxbtm6FZctg6dLy/0uXwhdfUADh2yW3F8O1o6ueDXfsCO3a+fbeUvn5cMUVSP03ctTouCfoyupFwgaUrCUmdFxFVmM1c6xs2QIFBWETc5kmTaBPH/jBD+Coo+jzyAMsWrGhartkr67wl7/ELlaRWlJvEraI1L2szIzw1cyZGbWTsEMTc2hyDk3M6enwve/BD38IffvCUUf5v+7dK3SuyujWNSHbJUWilSRD0iS+tLQ0BgwYQL9+/bjkkkvYvn17vEMiPz+fN998s9r18+fP5+677456f926daN///7079+fvn37csstt7Brl5/KeM2aNVx88cXVPvbbb7/l4Ycfjrj/k08+uSzuc889N+q4AObNm8eyZcvK7v/+97/n1Vdf3ad9SGR5uTn0692NtLQG9OvdjbzcHAqWF4WvZl5etO9PUFwM774LDz4Il1wCXbpAq1bw/e/D2LHw8MOwdq1PzFlZMG8efPYZbNsGixfDU0/Bb38L558PvXpV6Qk9ctRosqZMZ+KCrqTfbExc0JWsKdMTrtpTpDo6w64lTZs25f333wdg9OjRTJs2jRtvvLHGxxUXF9NwXy++j1J+fj7NmzcvS4SVn3fYsGEMGzZsn/a5cOFC2rVrx7Zt2xg/fjzjx49n5syZdOrUiWeffbbax5Um7GuvvbbKur1795KWlhbxx0VN5s2bx7nnnkvfvn0B+MMf/rDf+5Kqqqv67tS+DYtWbNy/y1+++w7++19YtMj//ec/PvmC7+B1yinQv3+1Z8z7IxHbJUWipTPsGPjBD35AYWEhzz//PCeeeCIDBw7kjDPOYN26dQBkZmYyfvx4zjzzTK644gpWrlzJD37wA4499liOPfbYssSVn5/PqaeeyqWXXkrv3r2ZPHkyOTk5DBo0iP79+7N8+XIANmzYwEUXXcQJJ5zAqaeeyhtvvMHKlSuZNm0a999/PwMGDOD1119nzJgx3HjjjQwdOpSbbrqJGTNmcP311wOwbt06hg8fzjHHHMMxxxxTY/Js3rw506ZNY968eWzatImVK1fSr5+flG3p0qUMGjSIAQMGcPTRR/PZZ58xefJkli9fzoABA5g0aRL5+fkMHTqUUaNG0b9//7J9ltqyZQvDhw+nb9++TJgwgZKSkirbPPvss4wZM4Y333yT+fPnM2nSJAYMGMDy5csZM2ZM2Q+IBQsWMHDgQPr378+1115bVivQrVs3br31Vo499lj69+/Pxx9/fGBvfD0WWvXdKK286rthA1+tvLAQ9uyFhYX+fkZmmGrm9eth7lz41a9g0CB/9nz66ZCZ6a9ZvvJKyMuDoiJYudLfjnDGLJJq6t8Z9i9+AcGZbq0ZMAAeeCCqTYuLi3nppZc466yzGDx4MG+99RZmxqOPPso999zDn/70JwDeffddFi1aRNOmTdm+fTuvvPIK6enpfPbZZ4wcOZJ33nkHgA8++ICCggLatGlDjx49GDduHG+//TYPPvggU6dO5YEHHuDnP/85v/zlLxk8eDDLli3joosuoqCggAkTJtC8eXN+/etfA5Cdnc2nn37Kq6++SlpaGjNmzCiL+4YbbuDUU09l7ty57N27l22lZzoRtGzZku7du/PZZ59x6KGHli2fNm0aP//5zxk9ejS7d+9m79693H333SxZsqSsFiI/P5+3336bJUuWhL1k6u2332bZsmV07dqVs846izlz5lRb5X7yySczbNgwzj333Crb7Ny5kzFjxrBgwQJ69+7NyJEjeeSRR/jFL34BQLt27Xjvvfd4+OGHue+++3j00UdrfN2pqGB5EYPHVVw2uDsUrd3EU08+Gf7yl/Xr4aWX4N//9mfQn37qH9ikCZx4Itx0EwweDCedFNfRo0SSRf1L2HGyY8cOBgwYAPgz7LFjx/LJJ59w2WWX8dVXX7F79+4KiWnYsGE0bdoU8AO/XH/99bz//vukpaXxaekXG3DCCSfQsaOfdbRnz56ceeaZAPTv35+FCxcC8Oqrr5a135aUlLBlyxa2bt0aNs5LLrmEtDBnKq+99hpPPPEE4NvjW7VqFdXrDjef+kknnURWVharV6/mwgsv5Igjjgj72EGDBlV7ffOgQYPo0aMHACNHjmTRokUR28ir88knn9C9e3d69+4NwKhRo3j88cfLEvaFF14IwHHHHcecOXP2ef/1VeVLtbp0qL7qu6ya2TnfKez55+Hkk+Gtt/yyNm189fbYsT5BH3ecT9oisk/qX8KO8ky4toW2YZeaOHEiN954I8OGDSM/P5/MzMyydQcddFDZ7fvvv59DDz2UDz74gJKSEtJDRkJqEvLF1qBBg7L7DRo0oLi4GPBJ+j//+Q9NmzatcWCQ0Oc9UFu3bmXlypX07t2bzZs3ly0fNWoUJ554Ii+88AI//vGPefTRR8uSb7SxVL6cqvR+6PJoRrYL94MiVGl5pqWllZVnqgvXXj0itxFjnmnMjEt2V+xhfe8f4LXXfJKePx8+/9zv5LjjfFX3eefBMcdUvJZZRPaLPkUxtHnzZg477DAAZs6cGXG7jh070qBBA5588kn2hht/OIIzzzyTv4RcR1r6w6FFixbVnmlXdvrpp/PII48AvhPYli1bIm6/bds2rr32Wi644AIOPvjgCus+//xzevTowQ033MCwYcP48MMP9ykW8FXiK1asoKSkhKeffprBgwcDcOihh1JQUEBJSQlz584t2766/X/ve99j5cqVFBYWAjBr1ixOPfXUqONIReHaq2eN2kNJoxblPayfb0dWj2MYOfEG3w79yCP+0qpp02D1anjnHfj972HgQCVrkVqiT1IMZWZmcskll/CDH/yAdu3aVbvdtddey8yZM/n+97/Pp59+us9nwX/+85955513OProoznhhBOYNm0aAOeddx5z584t63QWyYMPPsjChQvp378/xx13HEuXLg273dChQ+nXrx+DBg2iS5cu/PWvf62yzdNPP02/fv0YMGAAH3/8MVdccQVt27bllFNOoV+/fkyaNKnG13TSSScxefJk+vXrR/fu3Rk+fDgAd999N+eeey6nnXZaWVMBwIgRI7j33nsZOHBgWWc88GOBP/7441xyySX079+fBg0aMGHChBqfvz4JdzlWJNVdqrVm3SaW9O7H3gZpLPnya0YWFsKFF/qOZBs3wgsvwM9+BsGPVBGpXVZTlWGiOfLII90nn3xSYVlBQQF9+vSJU0SJRWOJR7Y/5ZPMx1eF6u0Kg4WEv/44Pz+f68ePYerpFWcqWlgIEx+DJe06wejRMHy47+mdYj238/PzGTJkSLzDSFgqn8jM7F3n3PH7+3idYYvUY9VdjpWVmVHtYzJu+BVjZzeqeKnWrDQyfnmTv+Tqnnt8z+4US9Yi8Vb/Op2JSJmCwlVhL8eqMhLZnj0wZw5s2cLIG34OzjExtwkFW3fRp+fhZD10lwYcEYkznWGL1FN5uTm0TDcWrai4vMJIZN9848+Ye/SAESP8nNC/+x0jCwtZsnkne0scSz4rUrIWSQA6wxapp7IyM7juZMfY2ZB9KWVt2KNzjT/9fiJcdx3MmAHbt8PQoX6s7hYtQG2QIglJZ9gi9Uhoj/BVRas4sj1knQ0T50L6ZP9/7VbHyF//Gh59FC691I8M+Npr/pppEUlYOsMWqSeqDngCY2bB3T+BJcGVdAsLYeLMBvCb38GECdChQ3yDFpGo6Qy7lpROr3nUUUdxzDHHMGXKlLIJK6qzcuVKcnNz6yhCqe/C9QifMQJ+94+Q3t7PNCbjkWw/CpmStUhSUcKuJaVDky5dupRXXnmFF198kdtuuy3iY5SwpTZVN+DJ5xuD6vAFXcj682OMvHJMXOITkQOTkgl7X0d+2lft27dn+vTp/OUvf8E5V+30mZMnT+b1119nwIAB3H///dVuJ6kt2uO1T88uYXuE9z2ii+/t/ekq9fYWSWIp14YdbmKDsTeOB6jVL7MePXpQUlLC+vXrad++fdjpM++++27uu+8+/u///g8g4jSbkpqiPl7XrCGjw+GMfWoV2T+l0qhmd8YpehGpTSmXsEPb+aB85KeJmRm1fvZROuxrpOkzQ0W7naSOGo/XnTthyhS4805G7tkDZ53HxFffp+Dz1RXnphaRpJdyCbtgeVF0Iz8doM8//5y0tDTat2/PbbfdVu30maEiTbMpqWlZ4SomboGC9dCnPWScARcfHRyvzz0Hv/41rFzpx/a+7z5G9ujByHgHLSIxkXIJ27fzVZzYoMLIT7Vgw4YNTJgwgeuvvx4zY/PmzXTu3JkGDRowc+bMsukzK08JWd12kprycnPo0MKYOtyVV3HPhqVroU+LJnDxxdC/PyxYAKedFu9wRSTGUq7TWUZmFmPnNqs4scHcZmRkZh3Qfnfs2FF2WdcZZ5zBmWeeya233gpUP33m0UcfTcOGDTnmmGO4//77D3iaTalfsjIzyBnlKk7ccSk89AZk7CyBBx+E995TshZJESl3hl3anjcxM4OC5UW11s4X6Wz4iCOO4MMPPyy7f9dddwHQqFEjFixYUGHbcNtJaipYXsTq46HfveVV4jedBlt2wcgVn0LXrvEOUUTqUMolbPBJWx1xJNF1at+GW17ayIwR5b2+x8yCToe2UbIWSUEpVyUukiwaUsKMEVQZuaxhA4t3aCISB/UmYZdeQiVSm+r6uCobJKWBsWrtN2FHLitau6lOYxKRxFAvEnZ6ejobN25U0pZa5Zxj48aNdXZ5XV5uDhm/vIapp69i593Qoy2R57IWkZRSL9qwO3fuzOrVq9mwYUO8Q4m7nTt36vrtCPa1fNLT0+ncuXMMIyqXlTGJ7At3lF1yePtZvs06tA3bj1x2YFc0iEhyqhcJu1GjRnTv3r3mDVNAfn4+AwcOjHcYCSthy+eJJyhY+VWFKvCRA6G4BIY9Dtv3mEYuE0lx9SJhiyStXbvgF7+AadPo07IJt/1zF/OWlF/GdUE/6NqlK0s+XRnvSEUkzpSwReLlq6/8kKL//S/85jcM3bKF7Jxp5I4urwIflQMXX35OvCMVkQSghC0SD++/D+edB99848cEv/BCFvbuRu5oKkz0kTsaJr7yYlxDFZHEUC96iYsklfnzYfBgf3vRIrjwQiCYmCbMZVy1PTGNiCQnJWyRuuIc3HcfXHAB9O0Lb78NAwaUre7SoU3Yy7i6dGhTp2GKSGJSwhapCyUl8POfw6RJfpat/Hzo2LHCJsUl/jKu0Ilpxszyy0VElLBFYm33bhg9mrypU+nXrgVpzz1LvwF9ycvNqbDZmvWbuONsmDgX0if7/3ec7ZeLiKjTmUgsbd0KF11E3iuvkNGmEdkXbw16gK9i7I3jgfIZ5Pr07ELnVqtYMqn84QsLNbKZiHg6wxaJlU2b4PTT4bXXyOrUluxL91Sc23r4drIyM8o2j9Vc7SJSP+gMWyQWvv4afvQjWLYM5syh4ILza+wBHqu52kWkflDCFqlt69f7M+vCQpg/n7yNX9My3Vi0wpVdYw3hJ/LQXO0iUh0lbJHa9NVXPlmvWgUvvACnnUZW725cd7Jj7GzIvrR8FLPRucafHlF1t4hERwlbpLZs2OCTdVERvPQS/PCHgK/2XjwOjurge36XjhO+dqvT2bSIRE2dzkRqwzff+DbrFSv8mfUPf0hebg79enfDOUefe/xmSybB3nth6nDo26trfGMWkaSiM2yRA7VlC5x1FhQU+GFHTz2VvNwcMm4cT/bw7Qwe56vASwdB6dxK81qLyL5TwhY5EDt2+Ek83nvPT+Lx4x8DkJWZQfbw7RUm8pgxws9t3bVLV/X+FpF9poQtsr/27oVRo+D11yE3l7xtW8nq3Y2C5UU45xg8ruLmg7vD9j2mua1FZL+oDVtkfzgH114L8+bBAw+QV7KXjBvHM/X0Vey8y9GjDWEn8tCoZSKyv3SGLbI/brsNpk+Hm2+GG24gq3e3ClXgt5/l26xnjCi/jEvt1iJyIJSwRfbVjBnk3XYbWW0OouCPd9Pn2VwKCldVqAIfOdB3MBv2uK8G16hlInKglLBF9kV+PnnjxpLRqgHZl35XNpHHhaurjmTWuZXvYKY2axGpDWrDFonWJ5/AhReSld6A7JElFSbyuO5kx+hc08QdIhIzStgi0di0CX7yE/KKi1n2XXGViTxuPdOPXDZxQVfSbzYmLuhK1pTpqgIXkVqjKnGRmhQXw4gR5K1aSUbrhvRo7DuRVZ7Io28vVX+LSOzoDFukJjfdBK+8Qlb71mRfvIvbz4Kxs1H1t4jUKZ1hi0Ty5JMwZQpMnEjBQ39hcHffbg3lE3k0awzTH1P1t4jEVkzPsM3sLDP7xMwKzWxymPWtzOx5M/vAzJaa2VWxjEckWnm5OfTr3om0K66gX8sm5J1wPH16dikbDGXkQD+Rx6s/8z3BlaxFJNZilrDNLA14CDgb6AuMNLO+lTa7DljmnDsGGAL8ycwaxyomkWjk5eaQ8ctrmHrWV+z8I0wdtYuMSf+PoT86h7Fzm6kqXETiIpZn2IOAQufc58653cAs4PxK2zighZkZ0BzYBBTHMCaRGmVl/pbsC3dUuGwre/h2Fr7yIllTpqsnuIjEhTnnYrNjs4uBs5xz44L7lwMnOueuD9mmBTAf+B7QArjMOfdCmH2NB8YDHHLIIcfNnj07JjHXB9u2baN58+bxDiNhRVM+p502lL7tfft0n/aQcQZcfDSkTzYWvPZaHUUaHzp+IlP5RKbyiWzo0KHvOueO39/Hx7LTmYVZVvnXwY+B94HTgJ7AK2b2unNuS4UHOTcdmA5w5JFHuiFDhtR6sPVFfn4+Kp/q1VQ+ebffTofmMHV4yBjgs2HpWujTq0u9L1sdP5GpfCJT+cRWLKvEVwOHh9zvDKyptM1VwBznFQIr8GfbInVv82ay7r6NnNFUrA6/FB5609RWLSJxFcuE/T/gCDPrHnQkG4Gv/g5VBJwOYGaHAkcCn8cwJpHwnIOf/YyC7XurjGI2uDts2enUVi0icRWzhO2cKwauB14GCoDZzrmlZjbBzCYEm90OnGxmHwELgJucc1/HKiaRaj35JHlPP03LptXMY92ra3ziEhEJxHTgFOfci8CLlZZNC7m9BjgzljGI1GjFCrj+erJaNuG6E3YxdravBi9twx6da/zpEVWHi0h8aaQzSW3FxXD55WBGwbbdLD4TjupQPopZn/Z+Ug9Vh4tIvGkscUlt994Lb7wBDz9Mlw5tWLSifBSzvff63uJdO7aNd5QiIkrYksI++ghuvRUuuQRGjaK4BMbMqjipx5hZUFwS70BFRFQlLqlqzx648ko4+GB4+GEwY836Tcy4rGJ1+B1nw5inN8U7WhERnWFLirrzTli8GP76V2jXjrzcHFo1bcCYWX71U6N8tXjnVtCnZ5f4xioigs6wJRV99BHccQeMGgUXXOAn+7hxPM9dvrfK6Ga5S5qRNUU9xEUk/pSwJfVcfbWvCn/wQQCyMjPIHr6dob386tLRzS56Mo2HpmtyDxFJDKoSl9Sybh288w5MnUreP1+mX+9uLPtsVdjRzTbvKFGyFpGEoYQtqePzz2HNGhg2jLziYjJuHM/U01fR99BqRjdT27WIJBAlbEkNzsG11/rbDz1E1m3l1eAZZ/g269DLucbObabJPkQkoagNW1LD00/Dyy/79uvOnSlYXsTq46Hfvf4Sri6t4Ypc+HIL9D2iK1lTslQdLiIJRQlb6r9vv4Vf/AKOPx7atwegU/s23PLSRmaMKB8zfMwsOKxDW5Z8ujKe0YqIhKUqcan/fvc72LDBX3MdaNgAZoyoOO/1jBF+uYhIItLXk9Rvixf7kcz+3/+DY48tW1z01cawPcOL1mpUMxFJTErYUn+VlMB110Hbtn6glEBebg4t0009w0UkqShhS/315JPwn//APfdA69Zli7MyM7juZFelZ/joXFPPcBFJWOp0JvXTli0weTKceCJccUWFVQXLi1g8TvNei0hyUcKW+umOO2DtWpg/n7xZeWRlZlCwvIjundvTqX0bFq3YyMiBfu5r8GfYExd0jW/MIiIRqEpc6p/PPoMHHoCrriLvs0/LRjTbeZfjb+eso8GerYzIbaSBUkQkqegMW+qf3/wGmjSBO+8k64ffrzKxx4xLdnP1/LZMXNCcguVF9OnZRQOliEjCU8KW+mXhQvLmzSPr0NYUHNaJZo0cqzdX3KT08q29e7+OT4wiIvtBCVvqj717ybtqDBktjezh31YYwaxhg/L2al2+JSLJSG3YUn/k5JD1VRHZo1yVEcx+9w+1V4tIctMZttQPO3bALbdQsJuwI5h9vgnSbza6d25P1pQ/qb1aRJKOzrClfvjzn8n74gtaNg0/t3XfXl3Zu7eER2fOUrIWkaSkhC3Jb+NGuOsuslo35bqTq85trRHMRKQ+UJW4JL8//hG2bKHAYPGZGsFMROonJWxJbl9+Sd4D95N1cDOaffcdi1agEcxEpF5SwpaklnfF5WQ0KSb70mJWb/aXcM0YQdklXWPnNiNriqrDRST5KWFL8iosJOuNhWRfRdlIZg0bwDXP+F7hfXt11QhmIlJvKGFL0sobM4ZVDs74q2+rzjjDV4VffLS/hGvJpyvjHaKISK1RwpaklHfvPWR88Abzrw6p/p7t13VooZHMRKT+UcKWpJR1121k/5QKk3pkX+qrw0saqt1aROofJWxJPh99RME326sd0SznqelqtxaRekcDp0jSyRs3jpbp1Y9opmQtIvWRErYklw8/JOv9t7nuFI1oJiKpRVXiklzuuIOC3RrRTERSjxK2JI+CAnj2Wbq0bsqiFTuqjGh29fy28Y1PRCSGVCUuySMrC5o1o7hxOmNmVawOHzMLikviHaCISOwoYUtyKCyEvDy49lrWfP0td5ztq8PTJ/v/d5wNa9ZvineUIiIxo4QtyeGPf4RGjcjr1o1WTRswZpZf/NQoWDIJOrfSYCkiUr+pDVsS3+rVMHMmeUNOJeMPk3ju8r0VRjdbuhZyl2iwFBGp35SwJfH96U9QUkJWYQHZw7dXGd3soifTeGi6BksRkfqtxipxM7vezA6ui2BEqti4kbxHHqJfq3SWrfgy7Ohmm3eUKFmLSL0XTRt2B+B/ZjbbzM4yM4t1UCKl8iaMJ6PRHqZe+h19Dw0/upnarkUkFdSYsJ1ztwBHANnAGOAzM7vTzHrGODZJddu3k/XCvLJJPjLOqDq62di5zTS6mYikhKh6iTvnHLA2+CsGDgaeNbN7YhibpLC83Bz69ezMsh0lTJwLeYv9IClZweVcTW6CiQu6kjVFbdcikhpq7HRmZjcAVwJfA48Ck5xze8ysAfAZ8JvYhiipJi83h4wbx5M9fHuVua5HDvTzXU9c0JUln66Ma5wiInUpml7i7YALnXOrQhc650rM7NzYhCWpLCszI2xv8IlzfbIeO1eXcIlI6ommSrx75WRtZk8COOcKYhKVpLSCwlVhe4MvW6dqcBFJXdGcYR8VesfM0oDjYhOOpLq83Jyyua5Lz7AhmOv6CFWDi0jqqvYM28xuNrOtwNFmtiX42wqsB/5eZxFKSsnKzOC6kzXXtYhIZdWeYTvn7gLuMrO7nHM312FMksIKlhexeJzmuhYRqazahG1m33POfQw8Y2bHVl7vnHsvppFJysnLzaFVEz8DV5/2/rrrkQP9GfbEBV3jHZ6ISFxFasP+FXAN8Kcw6xxwWkwikpSUl5tDxi+v4bkrnSb2EBEJI1KV+DXB/6F1F46kqqzMDLIv3KGJPUREqhGpSvzCSA90zs2p/XAkVRUUrmLwuIrLNLGHiEi5SFXi50VY5wAlbKkVkS7l0sQeIiJepCrxq+oyEEldoZdyZV9KWRv26FzjT4+o7VpEBCJXif/UOfeUmd0Ybr1zbkrswpJUkJebQ1ZmBss+WwVb4Cd9dCmXiEh1Ig1NelDwv0U1fzUK5s/+xMwKzWxyNdsMMbP3zWypmf1rH2KXJFY6wcfU01ex648wdTi8UOAv5dp7r7/ft5cu5RIRKRWpSvyvwf/b9mfHwRCmDwE/AlYD/zOz+c65ZSHbtAYeBs5yzhWZWfv9eS5JPprgQ0Rk39Q4+YeZ9TCz581sg5mtN7O/m1mPKPY9CCh0zn3unNsNzALOr7TNKGCOc64IwDm3fl9fgCSnguVFrN4M/e6FtEn+/+rNmuBDRKQ65pyLvIHZW/gz5bxg0QhgonPuxBoedzH+zHlccP9y4ETn3PUh2zwANMJPMNICeNA590SYfY0HxgMccsghx82ePTuqF5eKtm3bRvPmzeMdRo1GXHw+TfZuYcaI8k5mY2bBrrSWzHo2dkPVJ0v5xIvKJzKVT2Qqn8iGDh36rnPu+P19fDSzdZlz7smQ+0+Z2fXVbh3yuDDLKv86aIif+et0oCnwHzN7yzn3aYUHOTcdmA5w5JFHuiFDhkTx9KkpPz+fZCifpk0a8dgwKlSJzxgBV89vFNP4k6V84kXlE5nKJzKVT2xF6iXeJri5MOgwNgufcC8DXohi36uBw0PudwbWhNnma+fcd8B3ZvZv4BjgU6ReK1q7Keyc10VrN8UnIBGRBBfpDPtdfIIuPVP+Wcg6B9xew77/BxxhZt2BL/FV6aMqbfN34C9m1hBoDJwI3B9d6JLMuhx6MItWbNJAKSIiUYrUS7x7deui4ZwrDqrOXwbSgMecc0vNbEKwfppzrsDM/gF8CJQAjzrnlhzI80riy8vNYfuWbxkzi4pt2M805u4/q2e4iEg40bRhY2b9gL5AeumycJ3DKnPOvQi8WGnZtEr37wXujSYOqR+yMn/LrMtLWLu1fKCU7m2gpFEL9QwXEalGjQnbzG4FhuAT9ovA2cAioMaELRJOQeEXDB4HjdL8fNcAe/ZC+s1qvxYRqU6N12EDF+N7ca8Nxhc/BmgS06ikXuvSvBGLVlRctmgFdOnQJvwDREQkqoS9wzlXAhSbWUtgPRDNwCkiVa1bR/HO3YyZBQsL/Zn1wkJ/DXZxSbyDExFJXNEk7HeCIUT/hu85/h7wdiyDkvopLzeHfsd8jy/3wB1n+/br9Mn+/x1nw5r1qhIXEalOjW3Yzrlrg5vTgh7dLZ1zH8Y2LKlvSif7yB6+nYlzoXMrWDKpfP3CQl3SJSISSTRn2JjZhWY2BZgI9IxtSFIfhU72kXGGn/s6tEp87NxmZGTqki4RkepE00v8YaAX5WOJ/8zMznDOXRfTyCTplc53XbC8iGaNHKs3++WlPcMnzvWTffQ9oitZU7J0SZeISATRXId9KtDPBbOEmNlM4KOYRiVJL7QKfPC48sk9GjbwCXvkQD+N5sQFXVny6cp4hysikvCiqRL/BAhtXDwcPzKZSLVCq8AbpZVP7vG7f6gaXERkf0Sa/ON5/JjhrYACMyvtGT4IeLMOYpMkVrC8iMHjKi4b3B0+3wTpNxt9enZRNbiIyD6IVCV+X51FIfVOn55dWLRiVZXJPfr2UhW4iMj+qLZK3Dn3r9I/4GOgRfBXECwTqdbQH53DqJyKPcFH5fjlIiKy76LpJX4pfnKOfPxUm1PNbJJz7tkYxyZJbOErLzJ2UPnkHn3aw9hBMO+VF2t+sIiIVBFNL/EM4ATn3HoAMzsEeBVQwpYKQi/jcs6xeJwfwazUnr1w181F8QtQRCSJRdNLvEFpsg5sjPJxkkJKL+Oaevoqdt7l6NGGsBN8aDQzEZH9E03i/YeZvWxmY8xsDPAClea4Fql8GdftZ1Flgg9dxiUisv8iVombmQF/Bk4ABuPbsKc75+bWQWySRAqWF7H6eOh3b3mb9U/6wLDHYHuxLuMSETlQERO2c86Z2Tzn3HHAnDqKSZJQp/ZtuOWljcwY4a+3Lh3ZrHWDNLbuLY53eCIiSS+aKvG3zOyEmEciSa1hAz+SWeWRzRo2S493aCIi9UI0vcSHAhPMbCXwHb5a3Dnnjo5lYJJcir7ayODuFZcN7g5F326PT0AiIvVMNAn77Jo3kVSWl5tDy3Rj0QpXZWQz9QoXEakd1VaJm1l7M3sAeAiYAHzjnFtV+ldXAUriy8rM4LqTXZU5rkfnoF7hIiK1JNIZ9hPAu8BU4Fx8b/ExdRCTJJmC5UUsHgdHdag4stnabahXuIhILYmUsDs45zKC2y+b2Xt1EZAkn9KJPkrnuQZ/hj1xQdf4BiYiUo9E6iVuZnawmbUxszZAWqX7IoAm+hARqQuRzrBb4avELWRZ6Vm2A3rEKihJLproQ0Qk9qpN2M65bnUYhySxguVF3HR8xWVHtoeChZroQ0SktkRzWZdIRNWNctapvVpORERqi2bdkgNW7ShnOrpERGqNvlLlgBWt3RR+lLO1m+ITkIhIPRRVwjazwWZ2VXD7EDPrXtNjJHV06dAm7NzXXTqoSlxEpLbUmLDN7FbgJuDmYFEj4KlYBiXJpbjEVZn7eswsKC6Jd2QiIvVHNJ3OhgMDCS7pcs6tMbMWMY1Kksqadd8wY0TFy7ruOBvGPK0qcRGR2hJNwt4dzIvtAMzsoBjHJEmmS4vGdG61iyWTypctLNTEHyIitSmaNuzZZvZXoLWZXQO8CvwttmFJssh7Yibbd++qWiX+TGNN/CEiUotqPMN2zt1nZj8CtgBHAr93zr0S88gkKWRlTGLWGFi7tbxKvHsbKGnUQhN/iIjUohoTtpn9EnhGSVrCWfblhgpt10+NgouPhvSb1X4tIlKbomnDbomfrWsTMAt41jm3LrZhSTLIe/wxOjSHqcPLRzgbOxuWrlX7tYhIbauxDds5d5tz7ijgOqAT8C8zezXmkUnCy/rdTeSMrjjCWfal8NCbpvZrEZFati8jna0H1gIbgfaxCUeSRV5uDsu+/DrsCGdbdjq1X4uI1LJoBk75f2aWDywA2gHXOOeOjnVgkrjycnPI+OU19GhL2BHO+vTqGp/ARETqsWjasLsCv3DOvR/jWCRJZGVmkH3hDtZu9W3W2ZeGtGHPbUbWFFWHi4jUtmoTtpm1dM5tAe4J7lcYGNo5p27AKapgeRGrj4c/vgarvoGLZsK3O+GgxjD9semqDhcRiYFIZ9i5wLnAu4ADLGSdA3rEMC5JYJ3ateaWl76pMv91SXpbJWsRkRipNmE7584N/mtmLqmgYfFuHgvmv4by+a+vnh/fuERE6rNoOp0tiGaZpI6ib77T/NciInWs2oRtZulBu3U7MzvYzNoEf93w12NLKlq7li6Nq+kdrsFSRERiJtIZ9s/w7dffC/6X/v0deCj2oUkiyrt5MtsNTfYhIlLHIrVhPwg8aGYTnXNT6zAmSWBZc2Zpsg8RkTiIZrauqWbWD+gLpIcsfyKWgUkC+uorlm3Zpck+RETiIJrZum4FhuAT9ovA2cAiQAk7xeRNvkmTfYiIxEk0Y4lfDJwOrHXOXQUcAzSJaVSSkLLmzSbnp5rsQ0QkHqIZmnSHc67EzIrNrCV+EhANmpJq1qypUh2ecYavDtdkHyIisRdNwn7HzFoDf8P3Et8GvB3LoCTx5N08ufrqcE32ISISc9F0Ors2uDnNzP4BtHTOfRjbsCTRhFaHQ3l1+IUzjYcfVXW4iEisRZr849hI65xz78UmJEk4X31FwZZdmvtaRCSOIp1h/ynCOgecVsuxSILKu3kyLdN9NXjpGTZo7msRkboUaeCUoXUZiCSurLlPc90pVee+Hp1r/OkRVYeLiNSFaK7DviLccg2ckiLWrqVgyy4WnwlHdaBCL/G1W1UdLiJSV6LpJX5CyO10/DXZ76GBU1LD3Ln0aQy3/RPmLSlP1hf0A4pUHS4iUlei6SU+MfS+mbUCnoxZRJJYnn2Woc1akv32FnJHl1eHj8qBiy8/J97RiYikjGhGOqtsO3BENBua2Vlm9omZFZrZ5AjbnWBme83s4v2IR2Jlwwb4179Y2AhyR1cc4Sx3NCx85cV4RygikjKiacN+Ht8rHHyC7wvMjuJxafhpOH8ErAb+Z2bznXPLwmz3R+DlfQtdYu7vf4e9eynYuJXVm6HfveVV4jedBgXLi+IdoYhIyoimDfu+kNvFwCrn3OooHjcIKHTOfQ5gZrOA84FllbabCDxHxbZySQTPPgs9etBp+7fc8tImZoworxIfMws6tW8T7whFRFKGOedq3goIxhEvS/DOuYjzKQbV22c558YF9y8HTnTOXR+yzWFALv6a7mzg/5xzz4bZ13hgPMAhhxxy3OzZNZ7gp6xt27bRvHnzA9/R3r3wwQdw6KH89JfX8uRFWypcg72wEC5/riVPzf77gT9XHaq18qmnVD6RqXwiU/lENnTo0Hedc8fv7+OjqRIfD9wO7ABKAMNXkdc0AYiFWVb518EDwE3Oub1m4TYPHuTcdGA6wJFHHumGDBlSU9gpKz8/n1opn5kz4Ve/gv/+l6++3hJ2lLOvNm6tneeqQ7VWPvWUyicylU9kKp/YiqZKfBJwlHPu633c92rg8JD7nYE1lbY5HpgVJOt2wDlmVuycm7ePzyW17bnn4PDDyfvsM1qmG4tWuKqjnGkObBGROhNNL/Hl+J7h++p/wBFm1t3MGgMjgPmhGzjnujvnujnnugHPAtcqWSeArVvhn/+Eiy4i67YMrjvZMXa2rwbfs9f/H52rObBFROpSNGfYNwNvmtl/gV2lC51zN0R6kHOu2Myux/f+TgMec84tNbMJwfpp+x+2xNQLL8CuXXDRRRT8+UFuOh72FMOwx+G7XdD1YI1yJiJS16JJ2H8FXgM+wrdhR8059yLwYqVlYRO1c27MvuxbYui556BDBzj5ZDq1b8MtL23kiVEVe4gfdmjbeEcpIpJSoknYxc65G2MeiSSGHTvgxRfhiiugQQMaNoDHRlScB3vGCLh6fuTdiIhI7YqmDXuhmY03s45m1qb0L+aRSXy8/DJs3w4XXQRA0dpNYXuIF62NeFWfiIjUsmjOsEcF/28OWRbNZV2SjObMgTZt4NRTAd8TfNGKVeohLiISZzWeYQc9uSv/KVnXR7t3w/z5MGwYNGoEwNAfncOonIo9xEfl+OUiIlJ3NB+2lFu4EDZvhgsvLF/0youMHVRxHuyxg2CeJv4QEalTmg9bys2ZA82bw49+BEBebg6rilZxV6FP1E+NgpED/Zn2XTdr4g8Rkbqk+bDF27sX5s2Dn/wE0tPJy80h48bxzL+q/HKuscEQ7h1aqA1bRKSuRXOGXVnU82FLEnnjDVi/HoYPByArM4Ps4dsrXM6VfSlc8wyUNGxG1hSNciYiUpdiNh+2JJk5c6BJEzjHdyYrWF7E6uOrzoH9+SbIeWq6RjkTEaljsZwPW5KFcz5hn3kmtGgBUDbCWeU5sA87tK2StYhIHFR7WZeZ9TKzU5xz/wr5ewPobmY96zBGibV334UvvqjQO7xhAz+i2dBe0CitfISzhtEMtSMiIrUu0tfvA8DWMMt3BOukvpgzB9LS4LzzAN87/OuNGznjr75KPG+x30wjnImIxE+kKvFuzrkPKy90zr1jZt1iF5LUublzYcgQaNtWvcNFRBJUpISdHmFd09oOROKkoAA+/hgm+qv31DtcRCQxRaoS/5+ZXVN5oZmNBd6NXUhSp+bM8f8vuADwvcPDTfbx+SbImqLe4SIi8RLpDPsXwFwzG015gj4eaAwMj3FcUlfmzoXvfx86dQKqn+yjb6+uStYiInFU7Rm2c26dc+5k4DZgZfB3m3PuJOfc2roJT2Jq1SrfQzykd7gm+xARSUzRDE26EFhYB7FIXZs71/8fXl5hosk+REQS0/4MTSr1xdy50K8f9PL135rsQ0QkcSlhp6r16+H11+F3vwPQ5VwiIglOCTtVzZ/vhyQN2q91OZeISGLTQJOpas4c8tofQr9LhpGW1oBlhat0OZeISAJTwk5FmzeT98+XydjzDVNPL2LnXY4ebXw1eChdziUikjiUsFPRiy+SlVZC9mXFZZN73H6Wn40r9HKusXObkZGpqnARkUSgNuxUNHcuBbupUAU+ciAUl8Cwx2H7HqNPzy5kTcnS2bWISIJQwk41O3aQN//vtEz3Vd6hI5p1bgVdu3Rlyacr4xaeiIiEpyrxVPPKK2S53Vx3ir9sK7QKfHSuqQpcRCRB6Qw71QTV4YvPhKM6VBzRbO1WpypwEZEEpYSdSvbsgfnz6XPwQdz2z++Yt6Q8WV/QDyjqGu8IRUSkGqoSTyX//jds2sTQE39A9tswdTjsvNv/z35bE3yIiCQyJexUMncuNG3KwsKl5I6m7JKuob0gd7Sf+ENERBKTEnaqKCkhLyeHfk1gWeEXTJwLeYvLVw/uDgXLNcGHiEiiUht2isi7/XYyir8l+6dVJ/cYOdDf1wQfIiKJSwk7RWQ9NIXsn1Jlco+Jc/1sXGPnaoIPEZFEpirxei4vN4d+vbuybMOWsJN7LFsHExd01QQfIiIJTmfY9VjpHNfZw7czcUvVkc0WrYC+R2hkMxGRZKAz7HosdI7rjDOqjmymyT1ERJKHzrDrsYLlRQwe52+PHOj/T5zrq8H7HtFVk3uIiCQRJex6rFP7NvS5ZyMrNvnRzDLO8IOkTFyganARkWSjhF1P5eXm0GDPVv52SfllXGNmwc6SRjzwiKrBRUSSjdqw66mszAxmXLK7wmhmM0ZAs+YtVQ0uIpKElLDrobzcHFYVreKMv0K/e8tHNBvcHYrWbopvcCIisl9UJV7PLHj1VXIevZ/5V1Ud0axDC41mJiKSrJSw65m8Jx8tu5QLykc0u+YZKGmo0cxERJKVEnY9s+KL9WFHNPt8E+Q8pdHMRESSldqw65mO7VqwaEXFZYtWQNeObZWsRUSSmBJ2PVNc4i/fCh3RbMwsv1xERJKXqsTrmQ2btjBjhB/RrGC9HzDljrNhzNPqHS4iksyUsOuRvNwcWqZD51awZFL58oWF6h0uIpLsVCVej2RlZnDdyVUn+Rida5rkQ0QkyekMux4pWF7E4nFwVIeKVeJrtzp1OBMRSXI6w64H8nJz6Ne7G845+tzjly2ZBHvv9ZN99O3VNb4BiojIAdMZdpLLy80h48bxZA/fzuBx5ZN8FJf4tuyxczVYiohIfaCEneSyMjOqjGw2YwQMexy6dtGc1yIi9YWqxJNQaRV4WloDlhWuCjuy2fY9xpJPVypZi4jUE0rYSaa0Cnzq6avYeZejRxvCjmymy7hEROoXJewkE1oF3igNbj+r6shmY+c202VcIiL1jNqwk8yywlVM3FJ+yVbGGX4ks2GPw/Zd0L3zIWRNuV9V4SIi9YwSdhLJy82hQwtj6nBXYa7rUQOha9OmLDn0EPJnzmTIkCHxDlVERGqZEnYSycrMIGeUqzLX9YUzjYd374KrL45vgCIiEjMxbcM2s7PM7BMzKzSzyWHWjzazD4O/N83smFjGk+yWFa5i4lxImwT97oW8xb5H+JadjpF7S+Cyy+IdooiIxEjMzrDNLA14CPgRsBr4n5nNd84tC9lsBXCqc+4bMzsbmA6cGKuYkll11eFL10KflulwcAc44QT417/iHaqIiMRALM+wBwGFzrnPnXO7gVnA+aEbOOfedM59E9x9C+gcw3iSWmh1eKO08urwh940MrbtgksvBbN4hykiIjESyzbsw4AvQu6vJvLZ81jgpXArzGw8MB7gkEMOIT8/v5ZCTB4FhUUMHldxWWl1eMd77yO/Tx/Iz2fbtm0pWT7RUvlEpvKJTOUTmcontmKZsMOd7rmwG5oNxSfsweHWO+em46vLOfLII12q9YLOy82hVbMGLFqxt6zDGQQDpLRowpC//hU++QTMyM/PVy/xCFQ+kal8IlP5RKbyia1YVomvBg4Pud8ZWFN5IzM7GngUON85tzGG8SSl0pHNrv3+3irzXI+d05SMrbtgxAhVh4uI1HOxPMP+H3CEmXUHvgRGAKNCNzCzLsAc4HLn3KcxjCVphY5sFjrPdatmaTx04cWMfPJJn7BFRKRei1nCds4Vm9n1wMtAGvCYc26pmU0I1k8Dfg+0BR42f4ZY7Jw7PlYxJZu83Bw/uUfQdj1yoP/bsxfSby5h5GefwdFHQ9++8Q1URERiLqYDpzjnXgRerLRsWsjtccC4yo+T8qrw0sk9qrRdd+0Ib70Fd90VvyBFRKTOaPKPBFVaFX77WVRtu57bjIxBp/gNVR0uIpISNDRpgipY7i/japTm75e2XTdrDNMfm87Iu+6GU06Bbt3iGqeIiNQNnWEnqC4d2pTNcz1yICyZBK/+DNq1bcvIfv1hyRIYNSryTkREpN5Qwk5QxSVV57keM8svJycHGjb0o5uJiEhKUJV4glqzfhM/OxEumgnf7oDWTWHEAPjrfzdBbi6ceSa0axfvMEVEpI7oDDtBdWrfhhcK4LkrYdcf/f8XCqBT6+awejVcfnm8QxQRkTqkhJ2gGjaAGSOoMNnHjBHQsHg3tGgBw4bFO0QREalDStgJqmjtJgZ3r7hscHco2rILLroImjWLT2AiIhIXStgJqk/PLmW9xEstWgF9GgNXXBGXmEREJH6UsBPU0B+dw6icir3ER+XA0MYHwamnxjs8ERGpY+olnqAWvvIiYweVD5jSpz2MHQTzPkyDBvqdJSKSapSwE1RB4SoWj4M7zi5ftmcv3LVwa/yCEhGRuNGpWgLKy82hZbqFb8Pu2SU+QYmISFwpYSegrMwMrjvZVZn0Y3QOZGRmxTs8ERGJA1WJJ5i83BxWFa3irkLo0hquyIU1W30b9tptMHLU6HiHKCIicaCEnUBK58Cef5W/5nrRCj+15lOjoEMLmLiga7xDFBGROFHCTiClc2AP7eXvD+0F2ZfCNc9AiTUh60FVh4uIpCq1YSeQgsJVYUc3+3wjZD2QrepwEZEUpoSdICL1DO97aGtGjlayFhFJZUrYCSJyz/C74h2eiIjEmdqwE0TB8iIWj4OjOlQc3WztNhg5YUK8wxMRkThTwk4Qndq3oc89G1mxySfqsp7hL7SPd2giIpIAlLATQF5uDg32bOVvl5RfzjVmFuzcBQ/8bUq8wxMRkQSgNuwEkJWZwYxLdjO0FzRK85dzzRgBzZoepM5mIiICKGEnhILlRWEv5yr6dnt8AhIRkYSjhJ0AunRoE/Zyri4d2sQnIBERSThK2AmguMS3WYdezjVmll8uIiICSth1Ji83h369u5GW1oB+vbuRl5tTtm7N+k3ccba/nCt9sv9/x9l+uYiICKiXeJ0ondQje/h2Bo+DRStWMfbG8YCffatPzy50brWKJZPKH7OwUHNfi4hIOZ1h14HQST1Ke4FnD99OVmYG4Oe4Hju3WYUq8bFzm2nuaxERKaMz7DpQULiKweMqLhvc3fcOh2CO62XLmHjvnRTsgT69upI1JUuTfYiISBkl7Bgrn9TDlU2bCb4XeGiV98j3FjOydXtYtQrS0+MQqYiIJDJVicdYtZN65Fp5lff778NLL8ENNyhZi4hIWDrDjrFqJ/XY6sqrvO+6C1q2hOuui2+wIiKSsJSwYygvN4dWTRuQPnkvfdpDxhkwcqA/w564oKvf6JNP4JlnYPJkaN06rvGKiEjiUsKOkdJLuZ67fG/ZhB5jZ8PStZC7pBlZU4Lq8Dvu8NXgv/hFXOMVEZHEpoQdI6GXckFwKdelcNGTaTw0fbqvDi8ogNxc+NWvoL2m0RQRkeqp01mMFBSuCjuhx+YdJeVt13/4AzRtCpMmVd2BiIhICCXsGCi/lKvi8gqXcn3wATz9NEycCIccUvdBiohIUlHCjoGoLuW66Sbfyew3v4lrrCIikhzUhh0DBYWrIl/K9cor8PLLcN99cPDB8Q5XRESSgBJ2LQsd2WzkQH8ZF4RcyrV3rz+77tYNrr8+rrGKiEjyUJV4LauxOjw7GxYvhjvvhCZN4h2uiIgkCZ1h17KC5UXcdDzsKYZhj8N3u6DrwUF1+Jk/hiOPhFNPhREj4h2qiIgkEZ1h17JO7dtwy0vwxCjY9AdYMAFKHBx2aFv47W9h82Z46CEwi3eoIiKSRJSwa1nDBjBjBBXmvp4xAhqWFMOjj8LPfw5HHRXvMEVEJMkoYdeyorWbwg6YUrR+M3TvDrfdFp/AREQkqSlh16LyyT6g372Qt9gvX7QC+jQGHn8cmjePa4wiIpKc1OmslkSc7OO/kHXaj+GHP4x3mCIikqSUsGtJtZN9PA4Pte3EyOfmxDdAERFJaqoSrwV5uTksq26yj10wcsECaNYsPsGJiEi9oIR9gEqrwnu0IfxkH4e1g+99Lz7BiYhIvaEq8QNUWhW+dqtvs86+lPI27KcbkvWXB+IdooiI1ANK2AeoYHkRq4+HP74Gq76Bi2bCtzvhoMYw/dHHyue+FhEROQCqEj9ApSObTR0OO++G566Ew1tB64PbMPKnl8c7PBERqSeUsA9QtSObNdDQoyIiUnuUsA9QtSObrd0Un4BERKReUsKOUl5uDv16dyMtrQH9encjLzcHXnqJPk0sfO/wnl3iE6iIiNRL6nQWhdJLt7KHb2fwOFi0YhVjrx8D3xSTcfjhjH1uA9kX7SzvHT63GVlTsuIdtoiI1CNK2FEIO4rZZcVMfKY5Sz79FOY8x8TMDAqWF9GnZxeypmSpd7iIiNSqmFaJm9lZZvaJmRWa2eQw683M/hys/9DMjo1lPPurYHlR2Hbqgm++g/R0Ro4azZJPV7J3bwlLPl2pZC0iIrUuZgnbzNKAh4Czgb7ASDPrW2mzs4Ejgr/xwCOxiudA9OnZRe3UIiISV7E8wx4EFDrnPnfO7QZmAedX2uZ84AnnvQW0NrOOMYxpv2RkZjF2bjMWFsKevbCw0LdTZ2SqnVpEROpGLNuwDwO+CLm/Gjgxim0OA74K3cjMxuPPwAF2mdmS2g21ZgZthj3OYd/tovFBTdj93a7tX44a/dNNo0b/tK5DqUk74Ot4B5HAVD6RqXwiU/lEpvKJ7MgDeXAsE3a4kUPcfmyDc246MB3AzN5xzh1/4OHVTyqfyFQ+kal8IlP5RKbyiczM3jmQx8eySnw1cHjI/c7Amv3YRkREJOXFMmH/DzjCzLqbWWNgBDC/0jbzgSuC3uLfBzY7576qvCMREZFUF7MqcedcsZldD7wMpAGPOeeWmtmEYP004EXgHKAQ2A5cFcWup8co5PpC5ROZyicylU9kKp/IVD6RHVD5mHNVmoxFREQkwWgscRERkSSghC0iIpIEkiph1zTUaSows8fMbH3otehm1sbMXjGzz4L/B4esuzkor0/M7MfxibpumNnhZrbQzArMbKmZ/TxYrvIBzCzdzN42sw+C8rktWK7yCWFmaWa22Mz+L7iv8gmY2Uoz+8jM3i+9REnlU87MWpvZs2b2cfA9dFKtlo9zLin+8B3XlgM9gMbAB0DfeMcVh3L4IXAssCRk2T3A5OD2ZOCPwe2+QTk1AboH5ZcW79cQw7LpCBwb3G4BfBqUgcrHv14Dmge3GwH/Bb6v8qlSTjcCucD/BfdVPuVlsxJoV2mZyqe8LGYC44LbjYHWtVk+yXSGHc1Qp/Wec+7fwKZKi8/HHygE/y8IWT7LObfLObcC3xt/UF3EGQ/Oua+cc+8Ft7cCBfiR81Q+gPO2BXcbBX8OlU8ZM+sM/AR4NGSxyicylQ9gZi3xJ1TZAM653c65b6nF8kmmhF3dMKYCh7rg+vXgf/tgecqWmZl1AwbizyJVPoGguvd9YD3winNO5VPRA8BvgJKQZSqfcg74p5m9GwwZDSqfUj2ADcDjQZPKo2Z2ELVYPsmUsKMaxlQqSMkyM7PmwHPAL5xzWyJtGmZZvS4f59xe59wA/KiCg8ysX4TNU6p8zOxcYL1z7t1oHxJmWb0tn8Apzrlj8TMtXmdmP4ywbaqVT0N8c+UjzrmBwHf4KvDq7HP5JFPC1jCm1VtXOstZ8H99sDzlyszMGuGTdY5zbk6wWOVTSVBVlw+chcqn1CnAMDNbiW9yO83MnkLlU8Y5tyb4vx6Yi6/CVfl4q4HVQa0VwLP4BF5r5ZNMCTuaoU5T1XzgyuD2lcDfQ5aPMLMmZtYdP+/423GIr06YmeHbjwqcc1NCVql8ADM7xMxaB7ebAmcAH6PyAcA5d7NzrrNzrhv+++U159xPUfkAYGYHmVmL0tvAmcASVD4AOOfWAl+YWemMXKcDy6jN8ol3r7p97IF3Dr7n73IgI97xxKkM8vDTj+7B/0IbC7QFFgCfBf/bhGyfEZTXJ8DZ8Y4/xmUzGF+l9CHwfvB3jsqn7LUeDSwOymcJ8PtgucqnalkNobyXuMrHv9Ye+F7NHwBLS7+DVT4VymgA8E7wGZsHHFyb5aOhSUVERJJAMlWJi4iIpCwlbBERkSSghC0iIpIElLBFRESSgBK2iIhIElDCFokhM9sbzGy0xMyeMbNmCRBTppn9OszyTmb2bHB7SMhsVcMsmB3PzC4ws7778ZwP1DAqVuXtu5nZjmCIx4JglrErQ9afWzrbmEiqUMIWia0dzrkBzrl+wG5gQjQPMrOGsQ2rKufcGufcxWGWz3fO3R3cvQA/y1DUzKwN8H3nJ67ZF8udcwOdc33wA5n80syuCta9gB+VLO4/gETqihK2SN15HegVjBj1mJn9LziDPB/AzMYEZ+HP4ydYGGNm88zseTNbYWbXm9mNwWPeChIhZtbTzP4RTMjwupl9L1h+npn9N9j+VTM7NCSWY8zstWCO3muC7btZyDzrpYI4/mJmJwPDgHuDWoOeZvZeyHZHmFm4cbgvBv4Rst1KM7vTzP5jZu+Y2bFm9rKZLTezsD9onHOf46e9vCG47/BDq54bbeGLJDslbJE6EJwxnw18hB/d6DXn3AnAUHwCPCjY9CTgSufcacH9fsAo/JjNWcB25ycW+A9wRbDNdGCic+444NfAw8HyRfgz24H4sbF/ExLS0fhpJE8Cfm9mnWp6Dc65N/HDKU4Kag2WA5vNbECwyVXAjDAPPQWonMi/cM6dhP8RMwOf1L8P/CFCCO8B3wu5/w7wg5riFqkv6rzaTSTFNDU/nSX45JQNvImvzi1tR04HugS3X3HOhc53vtD5ub23mtlm4Plg+UfA0eZnJjsZeMYPpQ5Ak+B/Z+DpYMKBxsCKkP3+3Tm3A9hhZgvxPwjeZ989ClxlZjcClxF+Pt+O+GkHQ5XOA/AR0DzkNe4sHe88jMqzG60HavyhIVJfKGGLxNYO56ezLBNMUnKRc+6TSstPxE/JF2pXyO2SkPsl+M9vA+Dbys8RmApMcc7NN7MhQGbIuspjEu/vGMXPAbcCrwHvOuc2htlmB/5HSajQ11H5NVb3vTQQKAi5nx7sWyQlqEpcpO69DEwMEjdmNnB/d+T8fN8rzOySYF9mZscEq1sBXwa3r6z00PPNLN3M2uInuvhflE+5FWgR8vw78a/nEeDxah5TAPSKcv9hmVk34D78j5BSvfGTmIikBCVskbp3O9AI+DDo5HX7Ae5vNDDWzEpnUTo/WJ6Jryp/Hfi60mPexve0fgu43QXzHEdhFjAp6MjWM1iWgz9D/2c1j3kB/6NgX/UsvawLmA1Mdc6F/igYGuxbJCVoti4ROSBBW3wr59zvImyzCDjXOfdtLT3noUCuc+702tifSDJQwhaR/WZmc4GewGnOucpn8aHbnYhvz/+wlp73BGCPc+792tifSDJQwhYREUkCasMWERFJAkrYIiIiSUAJW0REJAkoYYuIiCQBJWwREZEk8P8Bl9Swc5X8sMEAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "cumul_prob = np.linspace(0.01,0.999,100)\n",
    "lognormal = stats.lognorm.ppf(cumul_prob,s = sigma, scale = math.exp(mu))\n",
    "\n",
    "plt.plot(lognormal,cumul_prob,color='red',label='Parametric Distribution',zorder=1)\n",
    "plt.xlabel('Permeability (mD)'); plt.ylabel('Cumulative Probability'); plt.title('Permeability Data Nonparametric and Lognormal Parametric Distributions')\n",
    "plt.ylim([0,1]); plt.xlim([0,600]); plt.grid()\n",
    "plt.scatter(df['Perm'].values,df[\"Perm\"].rank()/(len(df)+1),color='darkorange',edgecolor='black',zorder=10,label='Data')\n",
    "plt.legend(loc='upper left')\n",
    "\n",
    "plt.subplots_adjust(left=0.0, bottom=0.0, right=1.0, top=1.2, wspace=0.3, hspace=0.4); plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2cfa3584",
   "metadata": {},
   "source": [
    "#### Comments\n",
    "\n",
    "This is a simple distribution of heterogeneity metrics.\n",
    "\n",
    "Much more could be done, I have other demonstrations on the basics of working with DataFrames, ndarrays, univariate statistics, plotting data, declustering, data transformations, geostatistics, machine learning and many other workflows available at https://github.com/GeostatsGuy/PythonNumericalDemos and https://github.com/GeostatsGuy/GeostatsPy. \n",
    "  \n",
    "I hope this was helpful,\n",
    "\n",
    "*Michael*\n",
    "\n",
    "#### The Author:\n",
    "\n",
    "### Michael Pyrcz, Associate Professor, University of Texas at Austin \n",
    "*Novel Data Analytics, Geostatistics and Machine Learning Subsurface Solutions*\n",
    "\n",
    "With over 17 years of experience in subsurface consulting, research and development, Michael has returned to academia driven by his passion for teaching and enthusiasm for enhancing engineers' and geoscientists' impact in subsurface resource development. \n",
    "\n",
    "For more about Michael check out these links:\n",
    "\n",
    "#### [Twitter](https://twitter.com/geostatsguy) | [GitHub](https://github.com/GeostatsGuy) | [Website](http://michaelpyrcz.com) | [GoogleScholar](https://scholar.google.com/citations?user=QVZ20eQAAAAJ&hl=en&oi=ao) | [Book](https://www.amazon.com/Geostatistical-Reservoir-Modeling-Michael-Pyrcz/dp/0199731446) | [YouTube](https://www.youtube.com/channel/UCLqEr-xV-ceHdXXXrTId5ig)  | [LinkedIn](https://www.linkedin.com/in/michael-pyrcz-61a648a1)\n",
    "\n",
    "#### Want to Work Together?\n",
    "\n",
    "I hope this content is helpful to those that want to learn more about subsurface modeling, data analytics and machine learning. Students and working professionals are welcome to participate.\n",
    "\n",
    "* Want to invite me to visit your company for training, mentoring, project review, workflow design and / or consulting? I'd be happy to drop by and work with you! \n",
    "\n",
    "* Interested in partnering, supporting my graduate student research or my Subsurface Data Analytics and Machine Learning consortium (co-PIs including Profs. Foster, Torres-Verdin and van Oort)? My research combines data analytics, stochastic modeling and machine learning theory with practice to develop novel methods and workflows to add value. We are solving challenging subsurface problems!\n",
    "\n",
    "* I can be reached at mpyrcz@austin.utexas.edu.\n",
    "\n",
    "I'm always happy to discuss,\n",
    "\n",
    "*Michael*\n",
    "\n",
    "Michael Pyrcz, Ph.D., P.Eng. Associate Professor The Hildebrand Department of Petroleum and Geosystems Engineering, Bureau of Economic Geology, The Jackson School of Geosciences, The University of Texas at Austin\n",
    "\n",
    "#### More Resources Available at: [Twitter](https://twitter.com/geostatsguy) | [GitHub](https://github.com/GeostatsGuy) | [Website](http://michaelpyrcz.com) | [GoogleScholar](https://scholar.google.com/citations?user=QVZ20eQAAAAJ&hl=en&oi=ao) | [Book](https://www.amazon.com/Geostatistical-Reservoir-Modeling-Michael-Pyrcz/dp/0199731446) | [YouTube](https://www.youtube.com/channel/UCLqEr-xV-ceHdXXXrTId5ig)  | [LinkedIn](https://www.linkedin.com/in/michael-pyrcz-61a648a1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d2817dc7",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}
